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A novel expression for the near-contact pair correlation function of D-dimensional hard sphere 
systems is presented which arises from elementary free- volume arguments. Its derivative at contact 
agrees very well with our simulations for D = 2. For jammed states, the expression predicts that 
the number of exact contacts is equal to 2D, in agreement with established simulations. When 
the particles are wetted, they interact by the formation and rupture of liquid capillary bridges. 
Since formation and rupture events of capillary bonds are well separated in configuration space, the 
interaction is hysteretic with a characteristic energy loss Ea,. The pair correlation is strongly affected 
by this capillary interaction depending on the liquid-bond status of neighboring particles. A theory 
is derived for the nonequilibrium probability currents of the capillary interaction which determines 
j^jy the pair correlation function near contact. This finally yields an analytic expression for the equation 

of state, P = P(N/V,T), of wet granular matter for D = 2, valid in the complete density range 
. from gas to jamming. Driven wet granular matter exhibits a van-der-Waals-like unstable branch 

at granular temperatures T < T c corresponding to a first order segregation transition of clusters. 
For the realistic rupture length of the liquid bridge, s cr it = 0.07d, the critical point is located at 
T c = 0.274£ c b. While the critical temperature weakly depends on the rupture length, the critical 
density ci c is shown to scale with s cr i t according to s cr it = '4>j / '4>c — 1). The segregation 

transition is closely related to the precipitation of granular droplets reported for the free cooling of 
one-dimensional wet granular matter [J], and extends the effect to higher dimensional systems. Since 
the limiting case of sticky bonds, E c h 3> T, is of relevance for aggregation in general, simulations 
have been performed which show very good agreement with the theoretically predicted coordination 
K of capillary bonds as a function of the bond length s cr it. This result implies that particles that 
stick at the surface, s cr it = 0, form isostatic clusters. An extension of the theory in which the bridge 
coordination number K plays the role of a self-consistent mean-field is proposed. 
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PACS numbers: 47.57.Mg; 68.08.Bc; 83.80.Fg; 45.70.Vn 



I. INTRODUCTION 



' Dry sand trickles easily through chinks and crevices, as everyone knows well from the hour glass, or just personal 
experience. However, the addition of small amounts of liquid are sufficient to transform it into a plastic (or, more 
precisely, a viscoplastic) material. The same is true for all granular matter when a few volume percent of liquid are 
added, provided the latter wets the grains well and the grains are not too large. It is understood that this dramatic 
change, from a quasi-fluid to a solid behavior, is due to the formation of liquid bridges @, [H, 0, [E H S B S 03 G3 
f***. | between the granules wherever they come into contact. These liquid bridges mediate a cohesion force, and rupture as 
■ soon as the particle surfaces are separated by a distance s cr j t which scales as the cube root of the amount of added 
liquid [11]. These processes of formation and rupture of liquid bridges are the main cause of the observed dramatic 
changes in the mechanical properties of the material. Because of the generality of the effects, it has become common 
to study systems with spherical grains (usually glass beads), in order to ease theoretical modelling and to avoid side 
effects. We decided to follow this approach. 

In this work we show analytically that the peculiar interaction by capillary bridges gives rise to a first order 
transition, and we compute the critical density and the critical temperature. We shall focus on the two dimensional 
case, but many concepts carry over to dimensionality D = 3. Since there is no clear observation of a first order phase 
transition in the hard-sphere fluid for D < 2 [12|,|l3j], the added liquid leads to a qualitative change. More importantly, 
this transition is determined entirely by the geometric and energetic properties of the capillary bridges. 

A dry system of N hard spheres with diameter d confined to an area or volume V has no intrinsic energy scale, 
so that the equation of state is of the form P = T f(N/V) with the temperature T = (mviVi) and a nonlinear 
density dependence, /. The defined size of hard particles is conveniently used to restate the density n — N/V as the 



'Electronic address: axel.fingerle@ds.mpg.de 
^Electronic address: stephan.herminghaus@ds.mpg.de 



(A) 



cb 



-cb 



(B) 





-CD 



5 crit 



FIG. 1: The hysteretic interaction in wet granular matter. (A) Capillary bridges form at contact and mediate an attractive 
force _F c b. At the bridge length s cr it the bridge becomes unstable and pinches off. (B) This hysteretic interaction by capillary 
bridges gives rise to a well-defined loss of energy denoted by E c ^. The rupture length s cl -i t is largely exaggerated for illustration. 
While the particle diameter d is the only length scale for dry granulates, in wet granular matter there is a second scale set by 
s cr it. A realistic value is s cr it ~ 0.07d, which is realized when 1% of the jamming volume is added by a wetting liquid (with 
zero contact angle). Furthermore, the bond energy E c ^ defines an intrinsic energy scale, which is absent in dry granulates. As 
it is shown, the length and energy scale set by the capillary interaction give rise to a phase transition with a critical density <j) c 
and a critical granular temperature T c . 



dimensionless occupied fraction <f> = annd D /(2 D D) (<jd the surface of a Z?-dimensional unit sphere), which is the 
area fraction <p = jnd 2 for two, and volume fraction 4> = ^nd 3 for three dimensions. 

The capillary interaction of wet granular matter has a well-defined binding energy E c \^ (10|, and it has been 
demonstrated experimentally (T3 | under realistic dynamical conditions with impact velocities typical for strongly 
fluidized wet granular matter that the hysteretic character of the interaction is essential: the dominant mechanism 
of dissipation is the hysteretic formation and rupture of capillary bridges, the energy E c \, of which is irreversibly 
taken from the kinetic energy of the granular motion whenever a liquid bridge ruptures The bridge energy 

has been quantified [ll], [l4| , according to which £? cb is proportional dVW. Figure [T] illustrates this hysteresis in 
the Minimal Capillary model 11| applied here, which assumes a constant bridge force F^. This may appear as an 
oversimplification at first glance, but there is increasing experimental evidence that the details of the force law are 
insignificant for the collective dynamics on which we focus here (ill , [lil [l5| , as also confirmed from the point of view 
of dynamical systems theory [161 Il7| . 

Obviously, an external energy current has to be continuously injected to drive the system into a noncquilibrium 
steady state. In the equilibrium limit, E c h — > 0, we will have a pressure of the form P = T f{4>). It is the objective of 
this article to derive the equation of state for the hysteretic liquid bridge interaction of wet granular matter in such 
a driven state. In view of the intrinsic energy scale E'cb, this relation has to be of the form P = P(<fi, T/E c b). 

The equation of state is understood as an intrinsic property of homogeneous wet granular matter, kept in a stationary 
nonequilibrium state of granular temperature T. With this given temperature we may subsume various ways in which 
the system can be externally driven to compensate for the dissipation by rupturing liquid bridges, so that this granular 
temperature T is maintained over many particle diameters. 

We remark that in most experimental situations involving wet granular matter, the granular temperature is a 
nonlinear, even discontinuous, response depending on the details of the driving, such as boundary motion or air flow 
in air-fluidized beds. In this article we deliberately regard the granular temperature as the control parameter, so that 
the theoretical description of the boundary coupling is conveniently separated. Yet we emphasize that for the full 
description of an experimental situation one has to insert the equation of state into the equation for the external 
energy input, and then solve for the granular temperature as the nonlinear response to the external driving. 

We aim at describing the steady nonequilibrium states of wet granular matter, which are so multifaceted that at first 
glance one might think that aside from density and granular temperature further physical parameters are necessary 
in order to describe such a state. Yet as simulations have shown, states of wet granular matter far from equilibrium 
[l| are very well described by a single granular temperature T assuming a Gaussian velocity distribution, neglecting 
higher cumulants. Furthermore, it is known that the self-organized velocity distribution of free cooling wet granular 
matter has a vanishing fourth cumulant [181 ]. We point out that the condition of a locally isotropic and homogenous 
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state used in this work implies that the temperature field may vary only slowly over many particle diameters so that 
there is no strong influence by a heat current, which would otherwise be considered as a third parameter of the local 
nonequilibrium state. 

Throughout this study, we allow for a cert ain poly dispersity, < Ad/d < 0.1. (For higher polydispersity, the dense 
system undergoes a kinetic glass transition [ijj, [2(|). First of all, polydispersity is frequently used in simulations 
and experiments to prevent the monocrystalline state. Secondly, most systems of practical relevance exhibit some 
polydispersity. Another characteristic of 'real' granulates is that the surfaces of the grains are not ideal, bearing 
certain roughness. This does, however, only change the amount of liquid which must be added in order to achieve 
the capillary interaction: first some liquid is required to fill the crevices and tiny recesses in the grain surfaces, until 
the grains effectively have a smooth liquid coating, which is then completely wetted by all additional liquid. For glass 
beads, as those used in most of the experiments, this is typically the case above a volume fraction W m i n — 0.1% of 
liquid with respect to the jammed granular sample volume. We also require an upper limit on the volume fraction of 
the wetting liquid, so that the maximal length s cr i t of liquid bridges is of the order or below the polydispersity Ad 
of the spheres. This is to demand that s C iit/d w \/W /3 is smaller than Ad/d < 0.1, so that W < W m ax = 2.8%. 
This happens to closely coincide with the upper limit set on the liquid content to ensure that neighboring capillary 
bridge do not merge [2l| . For this range of the liquid content the capillary interaction is a truly pairwise interaction 
with the capillary force acting radially between pairs of particles. Another implication of roughness is that there is a 
substantial tangential friction between adjacent grains. This means that in principle one has to include all rotational 
degrees of freedom in the kinetic considerations for any statistical physical treatment of our system. However, we are 
here focussing on the effects due to the liquid capillary bridges, which mediate central forces. These do not couple 
to the (tangential) rotational modes. We therefore expect that the rotational degrees of freedom play, in our system, 
the role of a spectator heat bath which follows the translational dynamics, but does not influence it greatly, aside 
from a quantitative increase of the granular specific heat. In fact, recent experiments and simulations of wet granular 
systems [l5j show that this approach yields remarkable agreement with experimental data. In this work, we thus 
completely neglect all rotational degrees of freedom. 

II. DRY SPHERES AS STARTING POINT 

Before we add the wetting liquid to the hard sphere system, we investigate the dry case in this section and derive 
expressions for the pair correlation near contact, which will be extended to the wet case in the following section. 

Due to their finite size, the positions of hard spheres are not distributed independently from each other, as it is 
the case for the point-particles of the ideal gas. The configuration space of N spheres is not V , but restricted to a 
concave subset in which the systems moves chaotically as a high dimensional billiard. With the absence of an intrinsic 
energy scale, the dry system is athermal, which means that a change in temperature is equivalent to rescaling the 
time axis. The excluded volume gives rise to correlations in the particle positions, which are measured by the pair 
correlation function. Denoting by n = N/V the mean macroscopic particle density and by n m (r) = 5(r — r,) the 
microscopic density, the isotropic pair correlation G(r) is defined as the probability 

(M r ))particle at d Vo1 = n G (\ r \) d Vo1 = 11 d{s) d vol (1) 

to find the center of a particle in the shell d vol — <j D dr of radius r = ri+rj + s and thickness dr = ds centered 

around a reference particle. We have conveniently subtracted the particle radii + rj in the last equality of JT]), so 
that s > is the surface separation. The function g(s) is advantageous for polydispersity scattered around the mean 
diameter d [63| . because of its defined contact point, s — 0, which is smeared out in the function G(r). Furthermore it 
is the natural way to describe an interstitial liquid bridge between the considered pair of particles, with s the length 
of the bridge. For a certain liquid volume per particle and contact angle of the wetting liquid, there is a well defined 
critical bridge length s cr i t at which the bridge becomes unstable and ruptures. The mean density n is factored out in 
([T]) so that the dimensionless g would be equal to unity for all separations if there was no particle-particle correlation. 
Figure [3 shows the pair correlation of a fluidized state in which long range order is lost, so that g(s), respectively 
G(r), tends to unity for r ^> d. 

The forces in wet granular matter, hard-core repulsion and liquid bridge attraction, are short-ranged and radial, 
acting between pairs of particles over a separation range < s < s cr ; t with s cr it <C d. We are therefore interested in 
the short-range behavior of the pair correlation g(s) up to leading order in s/d. For such short particle separations 
the pair correlation g(s) is (up to a normalization constant) just the probability to find next neighbors at a separation 
s. Put in equivalent words: decomposing the pair correlation function g(s) — J^'kLi 9k(s) in contributions gu of the 
k's shell of Voronoi neighbors, we have g(s) = <?i(s) in the range of interest, < s < s cr ;t <C d. To shorten notation 
we suppress the subindex 1. 
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FIG. 2: The pair correlation of wet granular matter in a fluidized state resulting from a molecular dynamics-type simulation 
in D — 2 dimensions. The correlation function G(r) vanishes in the range (0, d) where the finite particle size leads to excluded 
volume. We use the function g(s) with the surface separation s of neighboring particles as it is convenient for wet granular 
matter where interstitial liquid bridges have the length s. Note that this is not exactly identical to the function G(d + s) shifted 
by one particle diameter d, since a realistic granular system has some polydispersity Ad around the mean diameter d. Aside 
from kinetic contributions, the pressure is due to the interaction forces which become dominant with increasing density. The 
internal forces in wet granular matter are short-ranged. Therefore our interest focuses on the sharp fall-off in the indicated 
range < s < s cr it of capillary interaction. This highlighted region indicates the typical range of s cr it, and corresponds to the 
region highlighted in Fig. 3] Furthermore, we derive more detailed correlation functions, g u (s) and g b (s), for unbound and 
capillary connected pairs, respectively, in order to describe the hysteretic interaction in wet granular matter. 



A. The Dense Limit 

Figure [3] gives an overview of results by [2(J, [22j (Fig. 15 therein), and [23[ for the phases of the two-dimensional 
system depending on density and polydispersity. For polydispersity below 0.1, there are two density regimes separated 
by the ordering transition at <j> [641 ]. These transitions are a purely geometric property (i.e. excluded volume effect) 
of the configuration space and are therefore athermal. To compute the radial next-neighbor distribution at densities 
above the critical density, <p > <j) , we consider the Voronoi tessellation of the system, which embeds each particle into 
a convex polygonal cell. The sizes {V^} of these Voronoi cells scales as (d + s) D , with the particle separation s. The 
mean cell size 

E™ K/N = V/N = 1/n is exactly the inverse density n. Hence, 



((-;)/ 



Ms -> 0) = </>(s ^ 0) (2) 



where the triangle brackets denote averaging over next neighbors which are in contact (s — > 0) with the center particle 
at the jamming density 4>(s — » 0) . We refer to those pairs of particles which come into contact at jamming as neighbors 
of type A, i.e the surface separation sa of A-neighbors is 

s A = at <j> = (f)j . (3) 

In the monodisperse limit for D = 2, 0j assumes the value of the triangular crystal, max = 7r/(2v3) = 0.91. 
Polydispersity decreases the (maximal random) jamming density <j)j and increases the critical density 4> a for the onset 
of triangular order as shown in Fig. [3] 



1. Contribution to the Contact Correlation: The A-Neighbors 

Since the Voronoi cells exchange their free volume, V — Vmin oc (l + 2)° — I7 an d the total volume is conserved we 
assume an exponential distribution of the free volume, which is well confirmed by experiments with dry granulates 
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FIG. 3: This plot reviews [20|, [23 ] (Fig. 15 therein), and [23 ] . The athermal transitions shown are properties of the configuration 
space of hard discs. Since the wet granular dynamics takes place in this configuration space, results for hard discs form the 
starting point for a theory of wet granular matter in two dimensions. At low (gas) and moderate (fluid) densities </>, the 
configuration space is probed ergodically and the system has low shear viscosity. As density is increased, the system gets 
trapped in a disordered state (glass for polydispersity above 0.1), or in a state with local triangular order. Both transitions, the 
glass transition (vertical line at c/> ~ 0.80) and the ordering transition (curved line ending at c/> = 0.71) can be detected by the 
rapid increase of the shear viscosity r\ (cf. [24j] for the ordering transition). While there is an athermal first order transition in 
three dimensions, it is at present discussed in the literature whether the transition region (double lines ending at <f> ) represents 
a fluid/solid coexistence (corresponding to a weak first order transition with a small jump of the entropy perparticle) or if 
there is an intermediate hexatic phase (according to the Kosterlitz-Thouless-Halperin-Nelson- Young scenario) p~3. . As the 
packing fraction <f> is increased further, the islands to which the system is confined in the configuration space shrink to points. 
This jamming limit can be detected by the divergence of the pressure p (at fixed granular temperature) under compression (for 
example using the particle expansion of the Lubachevsky-Stillinger algorithm). The maximal random jammed state (MRJ, for 
strict jamming as defined by Torquato, Truskett, and Debenedetti [251 ]) is the vertical curve at the right. Densities higher than 
MRJ are deep in the glassy regime. From the thermodynamic point of view, the system would cluster in phases separated 
according to the particle size, but this eutectic freezing-transition is kinetically suppressed [22| and unreachable. 



[26| . The conditions ^ and ([3]) determine the A-neighbor distribution uniquely: 
The volume element for D = 2 is 

dvol(s) = ob r {D - 1] dr = Trd (1 + s/d) As . (5) 

The contribution g& which A-neighbors give to the pair correlation is equal to the A-neighbor distribution Pa Q up 
to a prefactor, so that 

Ms)= 9 ? exp p^^ 1 ) (6) 

is determined as soon as we know the athermal contact value, g** = <?a(0). This contact value follows from the 
classical free volume theory (27j (which was based on [28|). 

in conjunction with the general relation between the particle- wall correlation gfjjaii and the pair correlation g^, 

^ = g^n = l + 2 D - V gf ■ (8) 
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As a consequence, we obtain 



,13-1 
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</>j/<t> - 1 



(9) 



close to jamming. Expression ((9]) is exact for D = 1. and has been confirmed as the asymptotic behavior of the 
diverging pressure close to jamming for D = 2 [24l [29j in event-driven simulation with accuracy 10~ 4 . We remark 
that this expression is not limited to weak polydispersity and has been confirmed for polydispersity far above 0.1 in 
the glass state [3(J HH . 

Inserting ([9]) in ([6]), we have as our first central result a closed expression for the near-contact pair correlation of 
neighbors which form exact contacts in the jamming limit (so-called A- neighbors) : 



9a(s) =5c'exp 



vD-l 



D 



4>& 



1 



- 1 



Eq. (JTUJ) implies for the derivative at contact, 

d 9a(0) 



'Vsi(o) 



(10) 



fill 



a quadratic dependence on the contact value = <?a(0). Eq. pip can be viewed as a consequence of normalization: 
the height of the contact peak is g** and so the width is of the order l/g^, which means that the negative slope is 
of the order g\. In fact, writing the A- neighbor correlation function 3a(s) in terms of the contact value g^, as we 



did in (|10p , is the natural form to express the density dependence of gA because this manifests that the coordination 
number of A- neighbors is density independent: 
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exp 
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- 1 



D-l 



ds 



(12) 



More significantly, Ka equals exactly the isostatic contact value 2D, which is obviously correct for particles on a line 
(D = 1) and is the accepted value for ideal discs and spheres in D = 2 and D — 3 dimensions respectively (22l l32l. [33l] . 
The finding (| 12(1 is an essential confirmation of consistency of our approach, since it is independent from conventional 
arguments based on the rank of the rigidity matrix (which accounts for global constraints on the degrees of freedom) 

As the contact value g^ © grows to infinity in the jamming limit, <f> — * <f>j, the constant integral (fT2|) implies that 
n 5a(s) becomes a delta distribution with 'weight' 2D at contact, s = 0. 



2. The Background Contribution: The B-Neighbors 

The configuration space is spanned by all particle positions {r^}. Consequently, a jammed configuration is - aside 
from a small fraction of rattlers [6f| - an isolated configuration point, and the set of jammed configuration is a set of 
discrete points. When the density is slightly relaxed, a finite system remains confined to a finite environment around 
the jamming point (cf. [33 |. p. 35). As density is lowered further, these environments are no longer isolated so that 
the system is able to migrate between theses 'islands of jamming'. 

The stability analysis of contact networks [13, HH, [H, [13, HH has put forth the result that frictionless spheres 
(except for the singular limiting case of a monodisperse crystal) jam strictly in an isostatic packing with 2D contacts 
per particle on average, as confirmed numerically [32| for D = 3, be the state random (glass regime in Fig. [3]) or 
locally ordered. Therefore we can identify within an island of jamming on average four neighboring particles in D = 2 
dimensions which are close to the reference particle, and which will be in contact with the reference particle, sa 0, 
in the jamming limit, <j> — ► 0j. These are the A-neighbors with the contribution gA to the pair correlation derived 
in (|10p . Furthermore, it is a mathematical fact that any discrete set of points in flat two-dimensional space has on 
average six Delaunay/Vorono'i neighbors [3t|, two of which have no contact to the reference particle, 5b (0) = 0. Hence, 
on the mean field level the following picture arises: Beside the four A-neighbors there are two B-neighbors which are 
sterically hindered by other particles from further approach to the reference particle. Summing up the contributions 
of A- and B-neighbors, 



g densc (s) = g A (s) + g B (s) , 



(13) 
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FIG. 4: The pair correlation near contact resulting from free-volume considerations. Close to jamming we distinguish between 
neighbors which form exact contacts in the jamming limit (contribution g\, Eq. (HQ}) and those that are blocked at positive 
separation s (curve gs, Eq. (HI}). The near-contact correlation is the sum of both contributions. For this plot the density is 
chosen to cj> — 0.8. The dashed curve sketches a typical second shell consisting of the second Voronoi neighbors. They are out 
of the interaction range, < s < s cr st, which is indicated by the highlighted stripe. 



gives us the pair correlation function near contact. The pair correlation near contact which arises from these blocked 
states, <7b, is discussed in detail in appendix [A] The essential result is that the configuration space of blocked states 
tends quadratically to zero in sb , so that to leading order the normalization of two B-neighbors for D — 2 determines 
the B-contribution in (II; 



Sb(sb) = WPb(sb) 
<j)cl \d) 



-[(l+f) 2 -!]/- 



1 + 



(t) 



(14) 



with cb = 4>max./4> ~ 1- I n Fig. [4] the resulting near-contact pair correlation 1131 for the dense regime is shown as the 
sum of <7a and <7b- 

B. The Dilute and Moderately Dense Regime 

In this part we turn to the free rheological regime, < <f> < <j) . When two spheres are closer than one diameter, 
s < d 7 they shield each other from certain collisions events. If one was to neglect three-particle correlations, the 
isotropic bombardment by 'thir d' p articles gives rise to the well-known attractive depletion force first proposed by 
S. Asakura and F. Oosawa Ef- As is evident from Fig. O summing up equal contributions over the accessible 
cross section is equivalent to the pressure exerted onto the submanifold indicated by the solid line in Fig. [5p and 
denoted by S. 

This depletion force, as well as the liquid bridge force which we will take into account in the next section, will 
affect the pair correlation function. A systematic way to study this effect has been worked out by Hansen et al. [12] , 
resulting in a Fokker-Planck equation for the two-particle distribution function. After integrating out the momenta 
and the center of mass coordinates, one finds that the depletion force as well as other non-entropic pair forces (such 
as the liquid bridge force), give rise to a Boltzmann factor, 



g(s) oc exp 



V[s] 
T 



(15) 



For the depletion force 



dcpl 



(16) 
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FIG. 5: Origin of the depletion force attracting neighboring particles that are separated by less than a particle diameter. One 
may either think of this as an entropic force, due to the decrease of excluded volume when the shells of excluded volume 
overlap. Equivalently one may view this as the net force due to isotropic bombardment. Obviously, the integration over the 
solid arc in A is up to a sign equivalent to the integration in plot B. In B the integration is over the outer solid arc, which is the 
configuration space of the third particle's coordinate at impact. Since the integration in B is projected by a cos- factor to give 
the axial symmetric force component, we can equivalently drop the cos- factor and integrate over the submanifold indicated by 
the solid line E in C. 



where 



n g c 



E ds 



conf 



-d lng(s) 



(17) 



is the infinitesimal logarithmic change of the excluded area (or the configuration space per particle, Vc 0n f), when the 
particles are separated by s < d, and E denotes the size of the corresponding section (line or area) in Fig. [5p. At 
contact, s = 0, the size of the integration section E is 



D - 1 



'V3 



(18) 



which yields Vdcpi to leading order in s. The depletion effect with the potential 



T 2 9 9c d 



s 

3d 



1 - tt-: - 



(19) 



for D = 3 has been confirmed in [43J by computer simulations. Polydispersity is known to have a minor effect on the 
depletion attraction [44J. For D = 1 (fT7)) gives the Poisson distribution Vd cp i/T = 4> flc* which is exact only for 
D = l. 

In two dimensions, the depletion potential is 



Vdepl 
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- 9 9c 

7T 



(4 arctan — ^ 
V W 

iff 



.1 + 

d 



7voiVF - C 



(20) 
(21) 



for D = 2 with 7 vo i(s) = 1 + s/d, the square root W(s) — y(T— s/d)(3 + s/d) and the constant C = 2tt/3 + y/3 to 
have Vdcpi = at s = 0. The first line (|2H)l is valid for < s < d, and the second line (|2"Tj) suffices for the region of 
interest, < s < s cr j t -C d. For the application of results on the near-contact decay of the pair correlation function, 
such as (|2"T1) . we prefer the exponential notation (used before in the dense case (fTU)) ) because it is most elegant to 
perform volume integration: 



gf exp 



(22) 



for D = 2. In this notation the dilute and dense behavior of the pair correlation are conveniently compared, showing 
that the result (|22|) for the gaseous/fluid regime differs by the factor l/4> max — 1.10 in the exponent from the dense 
result pi}|) close to jamming, so that according to (f2"2"]l the depletion force falls-off slower than the configuration density 
<j>g^ ■ We will now show that this is due to an over- estimation of the depletion force, caused by neglecting correlated 
three-particle events: when the plane of incidence of the third particle closely coincides with the symmetry plane E, 
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the incoming particle will hit in short sequence the pair of particles considered, which increases very effectively the 
exchange of momentum, i.e. the depletion attraction is reduced. 

To determine analytically and numerically the effect of correlated collisions which correct the Asakura-Oosawa 
result (|22[) we define the dimensionless measure 

4 fdcpi _ dg' c 

nnTdgf 5c 2 ' 1 ' 

for which the Asakura-Oosawa approach ([To]) and (|T5|) gives Zao = £ V3 » 2.205 (Line A in Fig. 0). Whe n we take 
correlated three-particle events into account, there are three contribution. Firstly, an attractive contribution Z\ > 
due to collisions on the front side of the pair, indicated by '1' in Fig. [51 which fall in the range — 7r/2 < p < ir/2. The 
corresponding value Zi is easily integrated. Isotropy of the state demands that the angle a between the symmetry 
axis of the pair P'P and the incoming momentum pi is uniformly distributed, as well as the impact parameter b (cf. 
Fig. [6j. These collision parameters are related by (a, b) = (tp + 9,ds'm9) to the position tp on P and the angle of 
incidence 9 with respect to the normal of P, which implies that ip is uniformly distributed and 9 is weighted by the 
cosine-factor cos 6*. Integration over —tt/2 < tp < ir/2 yields the axial force contribution 

Fa = 2Tng?d , (24) 

so that Z x = 8/tt ps 2.546. 

Secondly, the attraction is weakened by collisions hitting P in the remaining range n/2 < \tp\ < </5 max (s) (which we 
refer to as the 'broad side') giving rise to Z2 < 0. At contact p max (s = 0) is 27r/3. For these collisions the incidence 
is shadowed by the partner particle P' so that the angle of incidence 9 is restricted to — tt/2 < 9 < 9 msix (p). (Confer 
the collision event '2' in Fig. [6j) Some trigonometry determines 9 max (s,p) by the relation 

(1 +7voi (s) cos tp) sin 6> max = 1 - 7 vo i(s) cos6> max sin<^ , (25) 

which allows for an explicit function of ip at s = 0: 

2cos6 max (<p) = tan^ — y/l + 2 cos p . (26) 

After integrating over the impact momenta pi in the rest frame of P, the axial force imposed on P is 



F 2 {s) = -Tngfd \ dp cos tp / d0cos 2 0, (27) 

7T Jtt/2 J-tt/2 



where the cos</? projects the force on the symmetry axis of the pair PP'. The cos 9 factor appears quadratically in the 
integrand ([2T]) because of the cosine-distribution (or equivalent, because the Enskog collision frequency is proportional 
to the radial velocity (pi/m) cos 9) , and the transferred momentum which is pi cos 9 . Symmetry allows us to integrate 
over the upper half, ir/2 < tp < (p max (s) in (|2T[) and multiply by 2 with the general result 



f Vmax( S ) 

Zi{s) = ~z I dp cosp 

/tt/2 



7T 



1 max (s, f) sm29 max (s,ip) 

2 7T 2vr 



(28) 



and the numerical value ^(0) = —0.32813(9). 

Thirdly, the most obvious and important correction on the three- particle level comes from double collisions denoted 
by 3 in Fig. [SI The third particle hits first P' (gray arrow in Fig. [S]) from the broad side at tp' 6 (tt/2, 6* max ). The 
radial component pi cos9' of its incoming momentum p\ is transferred to P', which is why the third particle moves on 
tangentially to the circular cross section of P' with momentum p\ sin#' to collide shortly afterwards with particle P. 
Here the momentum transferred is the radial component with respect to P, p; sin^' cos#, so that 
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7T 



F 3 (s) = - Tngfd dp' cosp(p') cos6(s,tp') / d9' cos 9' sin 9' . (29) 



tt/2 ^0 



The collision point on P' described by <p'(<p) is related to p (the subsequent collision point on P) by cos(tp' — tp) = 
1 + 7voi(s) cos^j. The incident angle 9 on P is independent of 9' and given by sin#(s, tp') = I + 7 V ol(s) cost//. After 
the elementary ^'-integration we find 



'max(a) 

z *( s ) = Z¥ / dp' cosp(p') cos9(s,p') sin 2 9 ma .x(s,p') , (30) 

/tt/2 



Tt 2 
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FIG. 6: Three contributions to the effective force between a pair of particles P and P'. The collisions events 1 are attractive, 
while the events 2 cause a weaker repulsive forces. Furthermore the attraction is weakened by the temporally correlated 
collisions events 3. 



and ^3(0) = —0.091593(7). Summing up the three contributions gives Z COTT = Yli=i %i ~ 2.127 which is shown as the 
line B in Fig. [7] 

Based on our numerical data shown in Fig. we shall in the sequel assume the value 

Z s im — 2 . (31) 

By virtue of good statistics the simulation at <fr = 0.097 gave Z S i m — 2.0009 ± 0.0050, and Fig. [7] suggest this result to 
hold with few percent limits very well over the entire density regime < 4> < 4> considered in this subsection. The 
value Z = 2 determines the near-contact pair correlation uniquely to be 



g d ^{s)=gf exp (-cj>gf 



1+^ -1 



l + O 



(32) 



for D = 2. Satisfactorily, this result (f3"2")) has exactly the same functional dependence on the configuration density 
4> 5c' as the formula put forward for the dense case (fTUj) in the previous subsection. While three-particle collisions are 
obviously important since they shift Z in the right direction, no analytic explanation for this coincidence corresponding 
to the value Z = 2 is provided at present. Yet we shall see in the next section that any value other than Z = 2 would 
lead to inconsistencies when we introduce the liquid bridge interaction. 

We finally remark that the result (f3"2"| strongly differs from the 'Poissonian fluid' [45j, for which the contact 
correlation g** — 1 > is ignored. Even at the lowest density (<p = 0.1) considered in Fig. [7] the Poisson fluid would 
give Zp i sson ((/> = 0.1) = 1.7 which is 15% below the simulation value, and the deviation from Z = 2 grows with 
density <f>. 



III. THE PAIR CORRELATION UNDER THE HYSTERETIC INTERACTION 



In this section we dress up the pair correlation function in order to describe the status of the liquid-bonds which 
are created and ruptured hysteretically in wet granular matter. We will proceed in two steps: first, we introduce in 
part IHI Al the liquid bridges as hysteretic but forceless objects which follow the unperturbed particle dynamics. As a 
result, a direct relation of the dynamical system and the limiting case of isostatic granular packings [33, 38] at rest is 
found. I11 IIII Bl we turn on the liquid bridge force to its physical value, so that the bridges unfold their back-reaction 
on the granular dynamics. In the limit of low granular temperatures, T <C -E c b, the particles stick together. For this 
frozen state of wet granular matter the bridge coordination K is computed analytically as a function of the rupture 
length s cl -it, and we find very good agreement with simulations. 
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FIG. 7: Functional test of the near-contact pair correlation (|32l) . The vertical axis is proportional to the depletion force, 
fdepi oc — din g/ds at contact, divided by the configuration density cj> g^. This fraction Z is predicted to be density independent 
by (132[) and to assume the value Z = 2 (line C). Line A corresponds to the classical Asakura-Oosawa result, which is only 
valid for large beads immersed in a bath of small beads. In line B the corrections due to temporally correlated collisions events 
(derived in the text) have been taken into account. These events occur when a third particle of equal size strikes a pair of 
particles with a given separation s <C d as sketched in Fig. HJ] We proceed using the value Z — 2 (line C) because it agrees best 
with the simulation. Furthermore, Z = 2 corresponds to a near contact correlation function which is of exactly the same form 
as the function g\ we use in the dense case, when expressed in terms of the configuration density <f> gTi.^)- 



A. The Hysteretic Coupling 

Due to the hysteretic interaction, the pair correlation g is no longer a function of the particle separation s. In 
order to include the knowledge about the collision history the configuration space has to be enlarged in two respects: 
Obviously we distinguish between pairs with and without liquid bridges, which we denote by superscript indices, g h (s) 
and g n (s) respectively, for 'bridged' and 'unbridged' neighbors (cf. Fig. [5]). In addition, time reversal-symmetry is 
broken by the formation of the capillary bridge at contact. Hence we distinguish approaching pairs (with a negative 
relative velocity) which might collide and form a liquid bridge in the future, and those that move apart so that they 
can rupture the liquid-bond in the future. This relative velocity is denoted by a subscript arrow. 

As we discuss the radial pair distribution, contact and rupture become the important points on the s-axis of the 
pair correlation function. At theses points the functions g h and g u are coupled according to the hysteretic transition 
of the bond status. We use an intuitive notation to refer to these points: 

{The probability for a pair 
at rupture distance 
approaching without bridge. 




The probability for a pair 
at contact 

moving away with bridge. 



etc. 



The contact point s = (or to be more precise: the right-sided limit s = 0+) is denoted by a circled ©, and the 
rupture at s — s cr jt by the circled ©. Infinitesimally close to contact, there are four detailed correlation values: the 
bridge-connected and the unconnected states, either particularized by the sign of the relative velocity. The same is 
true for the left-sided limit s — s cr jt— of the rupture point. An infinitesimal distance beyond this point, at s — s cr i t +, 
there is only the unbound state possible with the two signs for incoming and outgoing velocities. This gives us in 
total ten detailed pair correlation coefficients. These are determined by the following ten equations describing the 
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hysteretic flow of probability, as it we can be read off from Fig. [8] 

Conditions on the contact shell 
9©^ = (33) 
5©->r = 9©^ + 9©^r (34) 

The Eq. (|33]) expresses that no particles rebound without a liquid-bond, but rather that all return with a bridge as 
stated by (|34l) . This implies that the particle number is conserved in collisions (in contrast to the absorbent dynamics 
modelled in for D — 1). 

Domain of capillary interaction 
9® = 7 b (s crit ) 9© (35) 
9® = 7 U (*crit) 9© (36) 

9c^© = 7pass 9c<—© 

(37) 

The functions 7 u (s) and 7 b (s) take into account the near-contact decay of the pair correlation without and with liquid 
bond, respectively. The last Eq. (|37|) describes spectator grains, i.e. grains which pass through the domain of possible 
capillary interaction without bridge formation. The fraction of these passing particles, 7 paS s = 1 — (1 + Scrit/d) 1 ^ 
(= 1/(1 + d/scrit) for D = 2) equals the gap between the considered cross section (2d + 2s CI it) D ~ 1 of the capillary 
interaction and the hard-core cross section (2d) 



Conditions on the rupture shell 

ScV© = (38) 

ScV© = .9c©^ (39) 

9^®+9c^® = .9"©^ (40) 

The Eqs. ([38l l39|) state that only unbound particles enter the domain of capillary interaction, and (|40|) describes the 
rupture of a capillary bridge when the pair escapes from the domain. 

The hysteretic capillary dynamics is coupled to the hard particle dynamics by the source term of new unbound 
pairs of particles entering the capillary interaction range: 



Source Term 

S©^r+ScV©/7 u (scrit) = (1-K/K sites )g$ (41) 

The left-hand side is the current of approaching unbound neighbors (measured at contact). If all neighbors were 
unconnected, K = 0, this current would equal the dry value g^.. But since there are K neighbors with bonds out 
of the -Ksites 'docking sites' which are sterically accessible for hquid bonds, the remaining unconnected fraction is 
l-K/K sitcs . 

The final tenth equation is the stationary state condition, which demands that the rupture frequency equals the 
binding frequency: 

Stationary state condition 

/bind = /rupt • (42) 

These frequencies follow from the probability to have a particle on the collision or rupture shell, respectively, mul- 
tiplied by the radial component of the relative velocity under the condition that the particle moves in the appropriate 
direction for the event to occur. This is analogous to the case D = 1 [l[, with the only difference that here we have 
to integrate over shells: 

./bind 
/rupt 

The volume factor 7 vo i(s) = (1 + s/d) 1 ^ ^ 1 ^ takes the increased size of the outer rupture shell as compared to the 
inner binding shell into account. 



i-D+l 



D 



T a u 

^ d 

T 

7T d 



and 



2 ° +lD \l- 3 7vol(Scrit) 



(43) 
(44) 
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FIG. 8: The hysteretic interaction in a wet granular gas or fluid can either lead to scattering or bound states. Note that the 
formation and rupture of the liquid bridge is spatially separated, which gives rise to a hysteretic loss and a coupling between 
the pair correlation functions g° for neighbors with and without, g u , capillary bridge. In this sketch the maximal liquid bridge 
length s C rit is drawn largely exaggerated. For a typical volume fraction of 1% wetting liquid added to the volume of jammed 
granular matter one finds s cr it/d ~ 0.07 [ill ]. 



Eliminating those correlation coefficients that are identically zero (|331 138|) , we can arrange the coupling equations 
for the domain of capillary interaction as a 6 x 6 matrix system: 



Collision: 


/ 


1 


1 


-1 











\ 


With Bridge: 





7 b 


7 b 





-1 





Unconnected: 




7 U 








-1 





-1 




Stationarity: 




-1 











7vol 







Spectators: 


V 











-1 





Tpass 




Source: 


1 














1/7 U 


/ 



S©<-r 
9©^ 

n u 

n h 

V s c V© J 



(1 - K/K sitcs ) 5 § 










(45) 



The 7-functions in the matrix are to be evaluated at s = s cr it- As it has to be on physical grounds, this system is 
non-singular with determinant (1 + s C rit/^) D 1 (2 + 7 pa ss) 7 b (scrit) > 0. 

The last row of the system (|45|) describes the creation of new liquid bridges as discussed before in the context 
of the equivalent Eq. (|41[) . We remark that here we used that the correlation <?jg of the dry system at contact has 
equal contributions from positive and negative relative velocities, immediately before and after the collision, which 
is still true for the wetted elastic particles we consider. This symmetry between positive and negative radial relative 
velocities is broken if one wishes to introduce a restitution coefficient < e < 1 to model inelastic collisions: the 
contact correlation of positive velocities is then increased by a factor 1/e as compared to the negatives. 

One should see clearly the very different meaning of K and -K s it es - The dynamical quantity K is the number of 
instantaneously existing capillary bonds: 



K = 2 D 



D 6 9 -f 
a 



t 

7 b (s) 7voi(s) di 



(46) 



K rapidly decays close to zero in dilute systems. As K comes closer to the value of i^sitcs in a very dense system, 
the binding frequency /bind oc g©<^ x oc K — -Ksites (02]) goes to zero because steric hindrance prohibits the formation 
of further capillary contacts: K — -Ksites gives the number of vacant sites for capillary bonds. Therefore -K s i tC s is the 
maximum number of 'docking sites' for capillary bonds. It is a pure geometric property and grows with s cr it, because 
Scrit > still allows for a slight rearrangement of particles in the formation of new capillary bridges without breaking 
existing ones. In the limit s cr ;t — > 0, -K s ites is the number of 'contact sites'. We therefore expect -K s it e s to equal the 
number of exact contacts, 2D = 4. So let us compute K s it e s{4>, Scrit) in the following paragraph. 

The maximum number of possible bonds, -K s it C s, is an athermal function of density (f> and the critical liquid bridge 
length s cr j t . We determine -K s i tC s from the obvious fact, that the granular dynamics is unaffected by the introduction 
of forceless bridges: for 7 b = 7" we recover the dry contact correlation + g^ = gj-j. This is the athermal limit, or 
high temperature limit of wet granular matter. 



High temperature limit 



9© + 9© 



„at 
9© 



(47) 
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FIG. 9: The capillary bridge coordination K in the low temperature limit, T <C -Ecb- As proven in the text, K converges 
to the athermal function Jfsites^t Scrit) in this low temperature limit. The solid line is K a ites(<j>o, Scrit) over a wide range of 
maximal bridge lengths s cr it . Points represent final states of free cooling simulations with 1000 particles of uniformly distributed 
polydispersity Ad = 0.06d. The open symbols are clusters with winding number one (cylindrical topology), connected over 
one periodic boundary on a rectangular domain. Such structures have internal tensile strength which necessitates a slightly 
increased coordination, visible as a small shift compared to the closed symbols which represent localized clusters (as the two 
examples drawn in the plot). As predicted by Eq. (|51[) of the presented theory, the structures emerging with exact contacts, 
Scrit — » 0, are found to be precisely isostatic, Twites = 4. The line is the the analytic result (|50|l . for which very good agreement 
is found with the simulations over the entire range of the capillary bridge regime, < s < 0.2r (with r the particle radius), 
which is indicated in the figure. Beyond this regime, the theory does not hold because in the derivation we limited ourselves to 
the leading order in s C rit / d. More importantly, the rupture length s C rit cannot be further increased beyond the capillary regime 
by simply increasing the liquid content in the granular sample. As mentioned in the introduction, liquid bridges residing on 
the same sphere would rather merge [2ll ] into more complicated objects. 
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FIG. 10: (A) The maximum number of 'docking sites' for capillary bonds on a particle, K s it cs , which is possible in two dimensions 
given the pinch-off length s CT it of capillary bridges and the density (j>. K s it as is independent of the temperature because it is a 
pure geometric quantity: the number of possible neighboring sites. Here -ffsites is shown as a function of <f> for different s cr it 
ranging from s cr it/d = 0.07 (solid curve), 0.04,0.02, to 0.01 (short dashes). As is shown in the text, A" s itcs is the coordination 
number in the zero temperature limit. While the mean coordination number K rapidly goes to zero with density for finite 
temperature, the zero-temperature limit is > 4 for all densities, because the system clusters. The dotted curve represents the 
limit Scrit — > 0. The longer s cr i t the closer .Ksites comes to six, the number of next neighbors in two dimensions. In the limit 
Scrit — » the coordination -Ksitcs converges to the number of exact contacts which is precisely four. (B) The capillary bridge 
coordination K drops down in the vicinity of the critical temperature as clustered structures break up. For this plot the mean 
density is chosen to be <f> = 0.75. 
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From the hysteretic bridge system (|45|) follows in this forceless or high granular temperature limit ([27] 

7 u 7v t 
1+7 P 



1 - i 7 , 7V °' ) ^sitcs = (1 + 7 u 7voi) K and (48) 



_at 

^ 1 + 7 U 7voi 

The 7-functions with the argument s suppressed are understood to be evaluated at s — s C rit- Inserting the Eqs. (|48[) 
and 09]) in gBJ) yields 

^ oDn , „at Jo"" 7 U (Scrit) 7vol(Scrit) ds/ri 

^ - 2 D ^ 9 ® l-7 u 7v ol /(l + 7p ass ) (50) 

- 4 1 - 7 - + 0(0" 4+0(Sc " t} ' (51) 

In the last line we have set D = 2, so that we could use (flQ|) and (fT4|) . The result ([ST) is the second important 
consistency test. Finding the number of exact contacts in the jamming limit to equal four in (|12p showed the 
consistency of the free- volume argument applied there. Here in (|5ip we find for any density that the different function 
.Ksites for the number of possible bridges sites equals four as well when s crit = 0. This is as intuitively expected and 
a confirmation of the consistency between the hysteretic system (|45j) and the near-contact pair correlation. In view 
of the numerical finding Z slm — 2 for the derivative at contact of the pair correlation (as defined in (|23|l ) we remark 
that the entirely analytic description by the hysteretic system gives in general K s n es = 8/Z + O (s cr it), which is why 
the consistency is non-trivial and the finding Z s - lm = 2 fits favorably into the entire picture. 

Thus the hysteretic system ([45)) provides a direct connection between the static granular properties captured in 
-Kgites and the granular system in motion at positive granular temperature which we are treating in general. We 
remind that K s n es is determined by the steric self-hindrance and therefore a pure geometric property independent of 
the granular temperature. When inspecting a snapshot of a close granular packing we can find local cases of contact 
coordination (s cr it = 0) higher than four. These are fluctuations within the granular ensemble, while K s it cs and K are 
mean-field quantities. Of course, for a finite bridge length s cr it > 0, a mean bridge coordination K S (0, if B ites) with 
-Suites higher than four is possible due to elongated bridges, as described by ([50[). Before we evaluated the expression 
dnO]) of i^sitcs f° r positive s cr i t (plotted in Fig. O, it is enlightening to switch on the capillary forces in the following 
section because this allows us to apply K s { tcs to 'frozen' wet granular matter. 



B. Switching On the Force of Capillary Bridges 

Under the attraction of a liquid bridge, the pair correlation g°(s) of connected neighbors falls off faster than g u (s) 
for unbound particles, depending on the granular temperature T / E c ^ compared to the bridge energy. The logarithmic 
derivative of the radial pair correlation is to be interpreted as the effective radial force O /3F = d s \ng(s), as 
discussed before in section III Bl This exponential dependence can be justified as the solution of the Fokker-Planck 
equation derived in [42]]. Moreover, in the context of the hysteretic interaction of wet granular matter this exponential 
factor has been successfully applied in the case D — 1 (cf. Eq. (3) in [TJ). Therefore we proceed by switching on the 
liquid bridge force to the physical value of the Minimal Capillary model Fb = -Bcb/scrit, including this exponential 
in the short-range dependence of the pair correlation function for bridges neighbors: 

7 b ( S ,T) = 7 u ( S ) expf-^ —) . (52) 

\ J Scrit / 

At low granular temperatures this exponential gives rise to shorter average bridge lengths, and describes the reduced 
probability that a bridge reaches its critical length s cr it. Therefore the hysteretic system (T4"5"]) describes the sticking 
of particles and the onset of clustering. 

We have discussed in the previous section IIII Al that steric effects in the dynamical system limit the mean number 
of bonds to a maximum of i^ s it eS ) and we derived that if sites converges to the number of isostatic contacts in the limit 
Scrit — ► 0. Here this connection is put on firm grounds with a clear physical interpretation attributed to -Ksitcs: -Ksitcs 
is the bridge coordination K of solid wet granular matter. 
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Proof of K — > K s ites in the low temperature limit. Solving (|45|) for K/K s i tos , we obtain 

K{T,s alt ,<p) _ 1 



(53) 



■^Csites (^crit 5 0) i + a(t)/t(t) ' 
with 

X(T) - 7 b (s crit ,T) K sitcs ( 7pass + 2) 7vol and (54) 
y(T) = 8/(T) cf>g$ ( 7pass + 1) , (55) 

where I{T) stands for the integral over bond states, 



Scrit 



TSr 



I(T) = / 7 b 7voi ds/d = / e-^TIT 7 u 7vol ds/d = — ^ + O (T 2 ) , (56) 
Js=o Js=o &cba 

which goes linearly to zero, while 7 b (T) cx e~ Ech / T vanishes for T — > faster than any power of T. Hence X/Y — > 
so that Eq. ([53"]) implies 

lim KT = KT site8 (57) 

as conjectured. 

This low temperature limit (T <C _E c b) is of general interest since it represents a sticky gas of ideal spheres, which 
serves as a model for aggregation in various areas of physics [48| and astrophysics 49]: once two particles had contact, 
the remaining degree of freedom is tangential motion. The analytic prediction of formula (|51[) is -Ksitcs = 4 in the 
limit of exact contacts, s cr it = 0. In order to evaluate ([50]) for positive s cr it we insert the near-contact decay 7 U given 
by the general results (fTT)]) . (|14[) . and ([3"2"]l , setting D = 2. The explicit expression for 7 U which we use throughout this 
article for results without free parameters is given in the appendix [B] Here we take into account known formulas for 
the contact value g@ at low densities, as well as higher corrections to the free volume theory. Inserting this expression 
in (|50p results in the curve shown in Fig. [51 We have performed simulations in this low-temperature limit. The wet 
granular matter was initially prepared in a gas state with T = 50-E CD and cooled by the formation and rupture of 
bonds. The insets in Fig. [5] show final states when the granular temperature T is more than one order of magnitude 
below _E c b and no further change in the configuration was observed on exponential time scales. The symbols in Fig. [9] 
have been measured in this final state. In perfect agreement with the prediction of Eq. (|?T|) . we find in the contact 
limit, s cr it — > 0, the coordination to be exactly 4. Moreover, the increase in the number of bonds per particle with 
the increase of the maximal bridge length s cr it is found to be in very good agreement with the simulations. 

Further analytic results for high densities are shown in Fig. [10] (A) . As is intuitively clear and shown by the family 
of curves in Fig. [TU] (A), the convergence of the limit s cr it — >■ is not uniform with respect to density, since -Ksitcs is 
pinned to the kissing number 6 of the monodisperse crystal density at ma x- 

IV. THE EQUATION OF STATE 

We are now in the position to derive the equation of state, P = P(T,<f>), for wet granular matter with capillary 
bonds tensile up to the rupture length s cr it- The cohesion of capillary bridges will reduce the pressure as compared 
to a dry hard-sphere system of equal temperature. By virtue of Eq. (|53p . we have the bridge coordination number K 
as a function of density <f> and granular temperature T . Since in the Minimal Capillary model [1 11 ] the bridge force is 
assumed to be independent of the bridge length s, the knowledge of the mean number of bridges K will allow us to 
evaluate the reduction of the pressure due to cohesion. Furthermore, the particle-particle collisions are enhanced by 
the bridge attraction, increasing the contact correlation. The contact correlation for wet granular matter derives 
from the Eqs. |@5J and (gBJ): 

wet at (1 + 7 pass)(l +7 b 7vol) ^ (c -o) 

J ® S ©(l+ 7 pass-7 U 7volK b +(2 + 7pass) 7 b 7vol/ U ^ ' 

with the integrals 

I u/b = f 7 u/b (s) 7voi(s) ds . (59) 
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FIG. 11: A local configuration of two-dimensional wet granular matter at moderate density. Since there is one particle in each 
Vorono'i cell, the mean area equals the inverse density. In section III Al we have used the Vorono'f tessellation to compute the 
derivative of the pair correlation at contact for a dry and dense system. For such a dense system, the Vorono'i cell resembles a 
hexagon with a size proportional to (d + s) 2 , where s is the particle separation. The cell borders are at half surface separation 
for polydisperse diameters (not half center distance), so that each cell contains one particle completely. 



The analytic expression (|58| for the contact correlation of wet granular matter, g©*, is indeed strictly greater than 
the one of the dry system, g^, to which it converges in the high temperature limit when the capillary energy E c ^ 
is small compared to the granular temperature T. This limit follows obviously from (|58[) because the functions with 
superscript index 'b' turn into those with 'u' for T E c b- In the low temperature limit, liquid bonds oscillate with an 
amplitude proportional to the kinetic energy which equals T on average, so that the probability to find the particles 
at contact, grows proportional to l/T, as can be derived easily from (|58]) using the expansion (j56|) . 



A. Frozen Degrees of Freedom 

As the system starts to cluster at temperatures close to £cb, voids remain between the clusters with linear dimensions 
large compared to the particle diameter. Clearly, this growing length scale, which is set by the sizes of clusters and 
voids, is not captured by the short-range behavior of the pair correlation function. Here we advance the theory beyond 
the level of two-particle correlations to take correlation on large scales, such as the collective particle motion in a 
cluster, in an approximative fashion into account. 

The collective motion of a cluster is due to stable capillary bonds which impose constraints, such that the internal 
degrees of freedom of clusters are frozen. Since K is the number of instantaneous capillary bridges of which the 

fraction erf E c ^/T^j with kinetic energies below E c \> forms stable bonds, we have 



frozen = K erf (60) 




for the number of frozen degrees of freedom. 

We are interested in the density of the remaining degrees of freedom. The idea is simple and powerful: As a general 
mathematical property of triangulations, there are on average precisely six Vorono'i neighbors independent of 
density or ordering. In Fig. [IT] we can observe that the Vorono'i neighbors with stable bonds contribute less to the 
area 1 jn of the Vorono'i cell. This picture suggest a two-fluid model with frozen and free neighborhoods as the two 
constituents. The fraction of frozen and free triangulation bonds is proportional to K{ Iozen and JsTf re e respectively, and 
the area contributions associated to each bond sum up to the total size of the Vorono'i cell: 

-^frozen + -f^frcc = 6 (61) 
rozcn . ^-free 6 / nrt\ 

+ = - . (62) 

^frozen ^frcc ^ 

Because of this reciprocal sum rule for the densities one may call this a reciprocal two-fluid model. 

The density rifrozen of the stable bond component follows analogously to ^ when averaged with the additional 
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exponential factor (I52|) due to the capillary force: 



/ frozen 



"'frozen 

/g" 1 * ... 7frozcn(s) 7vol(s) ds 
Jo 8 "" Tfrozen(s) 7vol(s) d 



7frozon(s) = exp 



/ „at 

-<P 9© 



(' + 5 



- 1 



Ecb _s_ 

-^crit 



(63) 
(64) 
(65) 



Without affecting the leading order in s/d one is free to replace the last s in the exponent ([65]) by s + s / (2d), so that 
the integral (|63p is elementary resulting in 



^frozen 

with a 



1 + 

^crit 

~d~ 



~d~ 

D 



1 1 

a e a - 1 
E c b d 



(66) 
(67) 



We point out that Eq. ((66)) implies Eq. (4) in [l[ for D = 1. 

From the Eqs. (|53"|) , (f6T))) - (|62[) . and (|66D - (IF71) follows the density of degrees of freedom which are not frozen out by 
capillary bonds, rif ree (T, s cr it , </>) • One may regard nf ree as the density of clusters. 

We remark that the two-fluid model of neighborhoods is the only concept presented in this theory of wet granular 
matter which cannot be generalized in a straight forward manner to three dimensions, because for D = 3 the number 
of Vorono'i neighbors (double counted per particle) is not a universal constant (such as 6 for D = 2 and 2 for D = 1), 
but depends on the granular order (reaching its minimum value 12 for close packing and its maximum of approximately 
15.5 in the ideal gas limit) The reason for this is that three-dimensional space cannot be filled with tetrahedrons, 
while flat space can be tiled by triangles. As a consequence, the number of constituents in the two-fluid model of 
neighborhoods would not be conserved for D = 3 and the numerator on the right-hand side of (|B"2"|) is not a constant. 



B. The Pressure of Wet Granular Matter 

Here we arrive at the pressure P(T,cj>) using the density nf rcc (T, <fr) (|61[) of degrees of freedom, the coordination 
K(T, <j)) (f5"3"]) , and the contact correlation g^ l (T : <fi) ([55)) . The pressure is the trace of the stress tensor 

P = -^trg. (68) 

The stress tensor g_ = <7 km + er force describes the flow of momentum. The kinetic term has components of™ = 
~ J2k (mv^v^ 5(t — r^)^. With the granular temperature T = (mviVi), its trace yields nT for uncorrelated 
particle motion (as in an ideal gas). In general we have the kinetic contribution 

P kin = n frcc T (69) 

wherein there frozen degrees of freedom have been taken out. For moderate densities, one may interpret (|69|) as the 
kinetic contribution to the pressure due to a gas of clusters. 

The interparticle forces F give rise to the Cauchy tensor a tmcc , which is the tensor product of the center-to-center 
vector r and the pair force F, 

gfor ce = ^free (p g ^ ^ (?()) 

so that <7 forco is diagonal for radial forces. The factor 1/2 assigns half of the momentum current to either of the 
interaction particles, i.e. r/2 may be seen as the transport vector within the Vorono'i cell. The Cauchy tensor ([70]) 
has contributions only by the unfrozen pairs of particles with density nf rGC , because in frozen neighborhoods the 
repulsive momentum exchanged in collisions is exactly balanced by the bridge attraction under the time average on 
the right-hand side of (|70p. 

A comment on the significance of the reciprocal two-fluid model as represented by Eq. (|61|) and (|62p is in order here. 
We consider for instance a compressed state of wet granular matter with K{ rozen around five and Kf ree around unity. 
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FIG. 12: The pressure P of wet granular matter is shown as function of the granular temperature T. The dimensionality 
is D — 2 and the covered area fraction is cf> — 0.1, so that at high temperatures the system is a dilute gas. The maximum 
bridge length is s cr i t = 0.07d. The behavior below the critical temperature T c = 0.274_E c b of wet granular can be understood in 
the following way: the system agglutinates to clusters. With these effective particles the pressure is reduced according to the 
reduced number density of effective particles. The breakup of clusters is reflected by the rising pressure around T c . The straight 
line is the athermal pressure of hard discs, P dly 

= n>9wa.llT which is reached asymptotically when the granular temperature is 
higher than the energy scale E c ^ set by the capillary interaction. 



While K{ ICC is small, the prefactor nf ree in (jTCT)) is not necessarily small. From (|61[) and (frj2")) follows that both, nf ree 
and nf r ozen, converge to nj as the system gets jammed (n — > nj), so that the repulsive dominated state is correctly 
described by the Cauchy tensor (|70|) which grows beyond all bounds as n — > nj. If one had (in contradiction to the 
additivity of areas) summed up densities linearly instead of the reciprocal sum rule (|62[) . the free density would vanish 
or could even become negative under such conditions. 

It is finally easy to determine the time average on the right-hand side of (|70[) for the two different forces acting in 
wet granular matter, the delta- force in collisions of hard particles and the flat force F c b = -E c b/s C rit of the capillary 
bonds. In a collision at time i co n the radial momentum Ap is transferred instantaneously: 



(F co u ®r) = (Ap®r S(t - t coU )) = 1 (Ap (r, -v) 9 ((r, -v)) S(r - d)) = -1 g^ t n<j D d D T 



(71) 



In the last equality the (^-function gives rise to the contact correlation and the trivial integration of angles leaves 
and 0-1 . 1 is the unity matrix and 8 is the Heaviside step function. Inserting (fTTj) in ([70)1 and taking the trace 
yields 



P, 



coll 



>D-1 



nft cc T g^ 



The cohesive virial due to capillary bridges is 



K 



E ch r (gi r 



Hence the final result 



P = n free T (1 + 2°- x g^) - n llcc E ch 



1 E ch 
n K — d 

L> Scrit 



K d 



2D Scrit 



(72) 



(73) 



(74) 



p.,., the contact correlation g^f and K have been derived 



where the last term is the bridge cohesion (|73[) . Since rif, 

explicitly in (frJT]) . (|55|) and ([55)1 as functions of <f> and T, we have the equation of state for wet granular matter, 
P = P(0,T). 

The Figs. [T2l and IT51 show the analytic result ([74")) as a function of the granular temperature T and the density <f>. In 
the high temperature limit wet granular matter behaves as a hard-spheres system. Below the critical point granular 
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FIG. 13: Isotherms of wet granular matter for the realistic rupture length s cr it = 0.07d. In the high temperature limit the 
liquid bridges forfeit their influence on the dynamics, so that the equation of state reduces to the hard sphere pressure. This 
can be seen by the two black isotherms of wet granular matter, of which the higher is at T = E c b and converges to the green 
curve in the limit T S> -E c b- The lower black isotherm is at T = 0.2i? c b and exhibits an unstable branch. The critical point is 
at T c w 0.274E c b (cf. Fig. [H for a close-up). 
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FIG. 14: A close-up of the transition region in wet granular matter. The dashed line in the main panel is the spinodal of the 
homogeneously driven wet granular system in D = 2 dimensions. The solid black lines are wet granular isotherms around the 
critical point, which is located at T cr i t = 0.273(5)i5 c b for s cr it = 0.07d. The change of the critical point with the amount of 
added liquid (represented by Scrit) is shown in Fig. 1151 The curve in the upper left corner is the athermal pressure P dry of the 
hard disc system [29t ] without liquid bridges, and the line at the bottom is the ideal gas pressure (P ld d D (4>) has a defined slope), 
pdry _ gat pid j g mcreasec ] compared to the ideal gas by the Enskog factor g^. The pressure of wet granular matter is reduced 
compared to the dry system P dry due to the capillary cohesion. The inset shows the spinodal in the temperature-density plane, 
where the critical temperature can be clearly determined. 
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FIG. 15: The influence of the rupture length s C rit on the position of the critical point in the phase diagram Fig. [14] of wet 
granular matter. The position of the critical point is described by the critical parameters (</> c ,T c ), which are plotted on the left 
and right vertical axis respectively. Solid lines result from the full theory (|74|) by solving for the intersection of d<f,P((fi, T) = 
and d^P(f>,T — 0. For the critical temperature we find a very mild variation with the rupture length, so that over the entire 
physically relevant range of capillary interaction we have T c ~ E c h/4. The influence of the rupture length s cr it on the critical 
density <f) c can be understood very clearly with the help of the dashed line. The critical density is such that the mean particle 

separation ~s = 4>j /4> c — l^j scales with the rupture length s cr i t . This shows that both intrinsic characteristics of the 

capillary interaction, the rupture length and the bridge energy E c h, determine the critical point of wet granular matter. 



clusters are predicted to segregate due to the mechanically unstable branch of the pressure as a function of density, 
which appears in Fig. [T2] below the critical temperature. Figure [H] provides a close-up of the critical point of wet 
granular matter and its spinodal. The critical density of this transition is high, because the particles have to be close 
enough in order to form a dynamical capillary network. As we show in Fig. 1151 the critical density is determined by 
the length scale of capillary bridges, such that the rupture length s cr j t scales with the mean particle separation s. 
Moreover, the rupture length is approximately four times the mean particle separation, s cr j t ~ 4s (dashed line shown 
in Fig. [15] for comparison). This result is to be compared with the very same ratio for the reported critical density 
of the un-clustering effect [l|: in the free cooling of dense one-dimension wet granular matter, the granular network 
was found to break up into granular droplets which precipitate out of the homogenous intial state, as soon as the 
density exceeded a critical value. This critical density was shown numerically and analytically to be set by s cr it ~ 3s 
[l|. The different prefactor is due to the additional cooling dynamics and the dimensionality D = 1. The theory of 
wet granular matter presented in this work predicts this transition to persist in higher dimensions. 

As we shorten the rupture length s cr it (which can be easily done experimentally by evaporating the wetting liquid) , 
the dry system is approached in such a way that the spinodal narrows in the T — <f> plane and is shifted to the 
jamming point, where it eventually shrinks to a line and vanishes. Figure [15] shows the convergence of the critical 
density to the jamming density. Since the capillary bridge regime sets an upper limit on the rupture length, the critical 
point is confined on the density axis between the ordering transition at 4> Q and the jamming density 4>j. The critical 
temperature almost exclusively depends on the bridge energy, according to T c w _E c b/4, over the entire capillary 
regime. 

With this discussion of transitions occurring in wet granular matter the presentation of our theory for wet granular 
matter is completed. The reader may find in appendix \C\ a brief methodical extension of the theory where a self- 
consistent equation is derived for future works. 



V. CONCLUSION 



Starting with the hard-sphere fluid, an expression (|10[) for the narrowing of the near contact pair correlation was 
derived, which describes in the jamming limit the delta-peak of 2D isostatic contacts per particle, in agreement with 
the excepted value of simulations. In the gas and fluid regime the fall-off predicted by this expression for the pair 
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correlation at contact was found to be well confirmed by simulations. We then addressed the nonequilibrium case 
of wet granular matter by the introduction of capillary bridges which are formed hysteretically. The description in 
terms of the pair-correlation function was extended with six different non- vanishing correlation coefficients which take 
the bridge status into accoount and allow for the hysteretic dissipative dynamics. The coordination number of bonds 
was computed analytically as a function of the rupture length of the capillary bridges, the granular temperature, 
and the density. The limiting case of strong bonds led to the sticky gas dynamics for which simulations have been 
performed which showed very good agreement with the analytic prediction of the coordination number. Based on 
the derived expressions for the contact correlation and the bridge coordination, we finally computed the pressure of 
wet granular matter analytically as a function of density and granular temperature. Here a method was put forward, 
which describes the effective degrees of freedoms in order to take the correlated motion of particles glued to clusters 
into account. The isotherms of wet granular matter were found to have an unstable branch which gives rise to the 
segregation of dense clusters. The critical temperature of this transition was derived to be approximately one quarter 
of the capillary bond energy. The critical density is directly related to the pinch-off distance of the capillary bridges. A 
close relation to the un-clustering effect reported in one dimension [l[ was shown, for which reason this effect persists 
also in higher dimensions. 

It will be interesting to probe the critical point of wet granular matter experimentally and by direct simulations. 
As we have shown, the position of the critical point is determined by the length and energy of the capillary bridges. 
These quantities can be controlled very accurately in an experiment of shaken wet granular matter. Therefore the 
measurement of the critical temperature will allow to discern between extensions such as the nonlinear coupling 
discussed in appendix [Cl 

Future analytic work includes the background contribution g# in the dense regime, since our numerics indicate that 
the pair correlation is flatter near the contact as predicted by <?a alone. This task might be addressed in conjunction 
with the analogous background contribution in three dimensions, for which in the jamming limit an integrable power- 
law divergence, gs oc l/s s , has been reported in numerical studies (with 8 = 0.5 51] or 5 — 0.6 [32| ) and experiments 
but is as well lacking a theoretical explanation at present. 

Acknowledgments Discussions with Martin Brinkmann, Svenja Hager, Jiirgen Vollmer, Klaus Roller and Mario 
Scheel are greatfully acknowledged. 



With <7a in (fit)]) we considered the four (cf. Eq. [T2"|) A-neighbors, which form isostatic contacts at jamming, sa — > 
for (f> — > ([3]), and are separated by sa according to Eq. ^ before jamming. Analogously, the separation sb of the 
two B-neighbors is weighted by 



While in Eq. ([?]) the denominator in the exponential is ca = 4>i/<j> — 1 so that sa — > at the jamming density, in 
Eq. (|A1[) the denominator is cb = </>max/0 — 1 since the blocked B is only forced to form a contact, sb — > 0, for a 
perfect crystal with <f> — > max . Of course this limit is kinematically unreachable because the system comes to rest 
at the jamming density 0j < max - '/'max would be reached. We note that cb is a small dimensionless quantity: for 
(f> > 4> = 0.71 we have < c B < 0.2774. 

Close to jamming, the B-neighbors are fixed in space by particles other than the reference particle. Except for 
arch-like constructions which are rare for frictionless particles, and would include second Vorono'i neighbors keeping 
B at a separation larger than our region of interest, sb > s cl -it, this hindrance is due to the A-neighbors. Therefore 
the probability <?b(sb) to find a B-neighbor at separation sb from the reference particle (sketched with hatching in 
Fig. 1 17[) is given by the integral over all configurations where four A-neighbors hinder two B-neighbors. 

The configurations will be weighted by a phase space factor C and the exponential factor Pb- We are above the 
ordering density <p a , so that the neighborhood has (by definition of the phase) hexagonal order as sketched in the 
inset of Fig. [TBI Projecting the configurations with the two B-neighbors blocked (gray subset in Fig. [T6j) on a single 
#-axis, we find the configuration space factor 



APPENDIX A: THE BACKGROUND CONTRIBUTION g B 



1. The Weighting Factors 
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FIG. 16: The angular configuration space of four neighbors close to the reference particle. These we denote as A-neighbors. 
The faceted inner subset shown in gray is the subspace conditioned to the property that two further particles, the B-neighbors, 
are hinder by the A-particles in approaching the reference particle. The projection of this subset onto an #-axis (for the 
angle between a blocking A-pair, 6\ or #3 in this example) gives rise to a linear configuration space factor C(8). Obviously 
a B-neighbor acts like a wedge driven between two A-neighbors, and therefore increases 6. This is taken into account by the 
weighting factor Pb(sb) which favors shorter separations sb between the particle B and the reference particle, depending on 
the density (f>. 



In the sequel we abbreviate 

7voi(s) oc 1 + s/d . (A3) 
for the volume factor (|A1[) in D — 2. Wide gaps of length sb are exponentially suppressed by Pb- 

2. The Configuration Space 

Let us now address the configuration space plotted in Fig. [T7] If the opening angle 9 of the A-neighbors exceeded 
0t(sa), 



COS^fA) = Vf^liA), (A4) 

2 a + s A 

the B-particle could slip through and turn into an A-neighbor, which is defined by having a free path towards the 
reference particle. This transition corresponds to the neck connecting different jamming island in the configuration 
space. Only along the line (PQ in Fig. [T7]) defined by 9 = 0c(sa), 

the B-neighbor can touch the reference particle, so that sb = 0. The Eqs. (|A4|) and (|A5jl define the upper boundary 
of the domain of integration for all sa, 

A >~\9 T ( SA ) s A /d>V2-l ' (A6) 



' max 



which is continuously differentiable but not smooth at the point Q. 
The lower boundary is 



9 s {s Al s B ) (s A + d) 2 + s| + 2ds B 

cos = — — — — , (A7) 

2 2(s A + d)(s B + d) { ' 
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FIG. 17: A section of the configuration space of neighboring particles. Within the gray domain the particle denoted by B is 
blocked: the two neighbors labelled A sterically hinder the particle B from approaching the reference particle (shaded). Only 
at the boundary 9c(sa) (curve PQ ranging from [sa,6]p = [0,2tt/3] to [s a ,0]q = [(V2- l)d,7r/2]) the B -neighbor can touch 
the reference particle. The probability <?b(sb) to find a B-neighbor at a separation sb follows from integrating over the gray 
domain, which grows with increasing sb- The lower bound, #s(sa, sb), is plotted for the values sb = 0.05, 0.30, and 0.60. Large 
areas spanned by this neighborhood are exponentially rare the higher the mean density cj>, so that the probability distribution 
in this plot concentrates in the vicinity of the upper left corner P as we come closer to the jamming limit. At the line QR the 
B-neighbor slips through and turns into an A-neighbor, so that QR is the transit to another jamming island in configuration 
space. The corresponding transition rate is proportional to the probability density along QR and therefore vanishes in the 
jamming limit. 




FIG. 18: The SA-#-plot of Fig. [17] with the full sb dependence shown on the additional vertical axis. 
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where B hits A. 

The simple lower bound on 6, 



i( s a) 



d/2 

l + S A 



(A8) 



which ensures that the A-neighbors do not overlap is without applicatory relevance, as it implies that the B-neighbor 
is pushed out to ss/d > \/3 — 1 ~ 0.73. This is suppressed in the dense regime 4> > 4> by the factor F of Eq. (|Al|) . 
The configuration space ends to its right in a cusp where the lower and upper bound intersect at 



2dsi 



2d 2 -d. 



(A9) 



This cusp converges to the point Q for sb — > 0. 

With the integration bounds l|A6p . (|A7[) . (|A"9|) , and the weighting factors (fA"Tj) . (|A2)) we have 



5b0b) = A/" Pb(s b ) 



A"P b (s b ) 
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e ma x(s A ) 
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d<9 C(0) 



(A10) 



(All) 



We emphasize that the configuration space (s A ,8) describes the relative position of one A- neighbor sketched sym- 
metrically in Fig. 1171 Since there are two independent A-neighbors involved, their configuration is the direct product 
(sai, $i) x (sa2i $2)- On this account the configuration integral is squared in (|A10[) . with the important consequence 
that the leading order in <7b(sb) is quadratical. The normalization constant M is determined by the knowledge that 
there are two B-neighbors. While the exponential prefactor dominates the long range decay, we expand the near- 
contact increase in s^/d. Substituting the dimensionless area z A = ((1 + SA/d) 2 — l) /ca for integration in favor of 
the particle separation sa, the expressions ij, i — 1, 2 are of the form 



h = ca 



1/CA 



v fi(c A z A ) dz A 



with 



h(x) 

h{x) 



3(x — l)a(x) 
2n 2 ^/(?>-x){x + l) 
2 2 6 VxTT 



6V3 



x-1 



y(x) a/3 — x a{x) v 'x + 1 



(A12) 

(A13) 
(A14) 



, . yjx + 1 
a(x) = 7r — 12 arcsin 



The integrals I{ can be treated by expanding the functions = f 



00 i/ CA 



(A15) 



=i/!-r(y+i,i/c A ) 



All incomplete Gamma functions can be eliminated by virtue of the recurrence relation (cf. (6.5.2) and (6.5.22) in 



I> + 1, l/c A ) = vT(v, l/c A ) + (-lyc^e- 1 '"* . 
As is apparent from the recurrence relation, the result will be of the form 

h =i? 4 (c A )+e- 1 / CA Si(c A ) . 
The regular part, for instance in first order of s^/d, 



Ri(c A ) = c A 
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FIG. 19: The radii of convergence r& for expansions around the jamming point. The contribution blocked B-neighbors give 
to the pair correlation can be expanded in a series around the jamming point, ca = 0. The radius of convergence is given by 
the Cauchy-Hadamard formula Tj = 1/ 4/ Kj for the term Kjis^/d) 1 . (A) The asymptotic divergence of the i?-series in (IA 181) 
poses no practical problem since few terms (less than 10) give sufficient accuracy. (B) The S-series in (|A18j) converges. 



is a series expansion about the point of jamming, ca = 0. It is asymptotically diverging due to the factorial which 
appears in the recurrence relation. Fortunately this does not restrain us from an excellent approximation, since for 
the relevant density, 4> > <f> Q , the quality of the expansion increases for more than 10 terms in the expansion (cf. panel 
(a) of the Fig.HU). 

The second part in (|A17|) . for which the first order of s^/d is given by 

a t x \/3 2 9 - 2V3tt 3 -27 + 2V3tt , A . 

Si c A ) = c A ^-+4 +4 — 2 + ... , A19 

27T 37T Z 37T^ 

and has a positive radius of convergence (cf. panel (b) in Fig. I19[) . This part is over-exponentially suppressed by the 
prefactor exp — 1/ca close to jamming. 

In the application to wet granular matter the sub-leading order (se/d) 3 < 4 • 10~ 4 is negligible for a realistic value 
of s < s cr it ~ 0.07<i, whereby we have the concise result 

9b(sb) = Afe~ ZB z B + O (4) (A20) 

with the abbreviation zb = ((1 + s-g/d) 2 — l) /cb and 1/cb ~ 05c' ■ The normalization 80^- J t/B dze = 2 according 
the two B-neighbors determines J\f in (|A20[) . Hence the result ([M]) . 



APPENDIX B: EXPLICIT EXPRESSIONS FOR THE PAIR CORRELATION IN TWO DIMENSIONS 

WITHOUT FREE PARAMETERS 



In our general derivation of the theory of wet granular matter we distinguished between the jamming density </)j 
and the (highest possible) crystalline packing max = 7r/(2\/3) achieved in monodisperse domains. The exact value of 
the jamming density </>j depends on many details such as the distribution of polydispersity and the jamming protocol 
for the increase of density. When we want to give explicit results without free parameters on the bridge coordination 
K(T, </>, s cr it) ^d the equation of state P = P(T, <fr, s cr it) we do this for weak polydispersity, where the difference 
between <f>j and 4> ma , x is negligible and the limiting case of 'dry' discs has been studied extensively. 



1. High Density 

For monodisperse 'dry' discs, <fij = </> ma x, there are higher order corrections to the free volume result ([9]) available 
in the literature which are incorporated in the final results on the bridge coordination and the equation of state for 
wet granular matter. These corrections are expansions with respect to x = </>j — <p fitted to simulations: 

5c donsG - (±+ao + a 2 .T 2 + ...) ^ = (± + a + a 2 x 2 + ...] fl + ± + ...\ , (Bl) 
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Scaled Particle Theory 
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TABLE I: The particle-particle correlation <?J* and the particle-wall correlation g^aii a ^ contact for different spatial dimensions 
valid up to moderate densities. The center column shows the results of the Scaled Particle Theory and the right column 
contains the exact exp ression for one dimension, and heuristic expressions [591 ] of Henderson [6fJ] for two dimensions and 
Carnahan-Starling [6l|] in three dimensions. 



Equation (|B1[) holds in the dense regime, <f> Q < cj) < i^ max , above (ft — 0.71. The numerical coefficients are ao = —1.07 
and 0,2 = 5.89 29], confirmed by our own simulations. Similar empirical expressions are also available for polydisperse 
discs in the glass state (Eq. (6) in [22]]). 



2. Low and Moderate Density 



For the analytic treatment an explicit expression for the contact correlation gf 1 in Eq. (|3"2")) is needed (as the 
counterpart to the dense expression (JbTJ) . Aside from the trivial one-dimensional case (66|, exact expressions for the 
contact correlation of hard spheres are unknown for the dilute regime. Yet there are well-established approximations in 
the literature resulting from Scaled Particle theory [&3.[55l , from the virial expansions [56] , as solutions of the Percus- 
Yevick closure [57], as well as heuristic expressions [581 ] such as the Carnahan-Starling formula with corrections to 
better fit simulation results (cf. Tab.QJ. 

As in the dense regime IB II we shall use the Henderson-Luding expression [2!| 



dilute 1 - 70/16 3 /128 



9c 



(1 



(1- 



(B2) 



for the uncaged regime, < <f> < 4> 0l and the merging function m(cf>) = 1/(1 + exp((0 o — (f>)/mo)) with a cross-over 
width mo = 0.0111 to smoothly connect the dense (1B1[) and dilute (|H2|) expressions [29| : 

(B3) 



g(s) = m(« ff dilutc ( S ) + (1 - m{<f)) 5 dcnse (s) 



with ,g dilutc (0) = gf lutc and 5 dcnsc (0) = g dcnso as given by the Eqs. flBT]) and p2| . The near-contact decay has been 

(B4) 
(B5) 



established in the Eqs. (jlO]) and dT4J) for (j>> <j) , and in Eq. (J32J) for < < O : 

<Ulutc (*) = 9? 7dilute ^d 



9 



dense 



(s) 



9c Tdense 9 c Tdilutc 



i 



up to leading order in s crit with 



7 U ( S ) 



dilute 



exp 



.9c 



(B6) 



With the contact expressions (|B1I IB2[I . as well as the short-range decay formulas (fTUl [3^|) . we have sufficient 
information on the dry system over the entire density range. We may therefore proceed by introducing the hysteretic 
capillary bridges. 
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FIG. 20: A typical graphical solution of the self-consistent equation (|C2[l . Here the density is chosen to be 4> = 0.6 and the 
granular temperature is T = 0.2_E c b. 



APPENDIX C: SELF-CONSISTENCY OF BRIDGE COORDINATION K 

All results presented so far on the coordination K((j>, T) and pressure P{4>, T) allowed explicit analytic results. Here 
we want to demonstrate how to treat more complicated source terms of the hysteretic system (|45|) numerically. Such 
an extension of the theory could be motivated as follows. The current of free (unbound) approaching particles could 
be a function of the free density Uf ree instead of the mean density, since some of the unconnected neighbors traverse 
the voids between clusters, so that Eq. (|4T|) is changed to 

cj) g^_ t + $ ,g c V©/7 u (scrit) = (1 - K/K sitcs ) froc g$(cb hcc ) . (CI) 

Obviously this approach is a lower estimate for the current of freely approaching particles, which is why (|C1[) is 
considered as a methodical example rather than a physical competitor to the theory presented above. 

With the altered Eq. (|Clj) the hysteretic system (j45|) can still be solved analytically to find the correlation coefficients 
g = {5© 4 - r >5© 4 _r>fl'©-^r)5c-K£))fl , c-^©)5S-©}- Unlike before, due to the coupling JUT]) and the Eqs. ((60l) - d62j) , the 
correlations g are a highly nonlinear function of K . Therefore Eq. (|46[) becomes a nonlinear self-consistent equation: 

if(g(/C,0,T),0,T)=/C. (C2) 

The physical value K(4>, T) of the coordination is the solution JC of (|C2j) . The numerical solution of (|C2j) is found to be 
very robust, as Fig. l20l indicates . Plugging the resulting self-consistent K(4>, T) back into the equation for the pressure 
(f7"4")) of wet granular matter, we find that the critical point is shifted from T c = 0.273(5)i? c b to T c = 0.216(5)£ , c b- 
This reduction of the critical temperature is intuitively clear since with less particles arriving to form bonds, the wet 
granular matter 'evaporates' at lower granular temperatures. 
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